function result = mie_ab(m,x)

  % Computes a matrix of Mie Coefficients, an, bn,
  % of orders n=1 to nmax, for given complex refractive-index
  % ratio m=m'+im" and size parameter x=k0*a where k0= wave number in ambient
  % medium for spheres of radius a;
  % Eq. (4.88) of Bohren and Huffman (1983), BEWI:TDD122
  % using the recurrence relation (4.89) for Dn on p. 127 and
  % starting conditions as described in Appendix A.
  % C. Mätzler, July 2002

  z = m .* x
  nmax = round(2 + x + 4 * x.^(1/3))
  nmx = round(max(nmax, abs(z)) + 16)
  n = (1:nmax)
  nu = n + 0.5

  sx = sqrt(0.5 * pi * x)
  px = sx .* besselj(nu, x)
  p1x = [sin(x), px(1:nmax - 1)]
  chx = -sx .* bessely(nu, x)
  ch1x = [cos(x), chx(1:nmax - 1)]
  gsx = px - i * chx
  gs1x = p1x - i * ch1x
  dnx(nmx) = 0+0i
  for j = nmx:-1:2      % Computation of Dn(z) according to (4.89) of B+H (1983)
    dnx(j - 1) = j ./ z - 1 / (dnx(j) + j ./ z)
  end
  dn = dnx(n)          % Dn(z), n=1 to nmax
  da = dn ./ m + n ./ x
  db = m .* dn + n ./ x

  an = (da .* px - p1x) ./ (da .* gsx - gs1x)
  bn = (db .* px - p1x) ./ (db .* gsx - gs1x)

  result = [an; bn]
